Dielectric Response of Periodic Systems from Quantum Monte Carlo Calculations 
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We present a novel approach that allows to calculate the dielectric response of periodic systems 
in the quantum Monte Carlo formalism. We employ a many-body generalization for the electric 
enthalpy functional, where the coupling with the field is expressed via the Berry-phase formu- 
lation for the macroscopic polarization. A self-consistent local Hamiltonian then determines the 
ground-state wavefunction, allowing for accurate diffusion quantum Monte Carlo calculations where 
the polarization's fixed point is estimated from the average on an iterative sequence, sampled via 
forward-walking. This approach has been validated for the case of an isolated hydrogen atom, and 
then applied to a periodic system, to calculate the dielectric susceptibility of molecular-hydrogen 
chains. The results found are in excellent agreement with the best estimates obtained from the 
extrapolation of quantum-chemistry calculations. 

PACS numbers: 77.22.-d, 71.15.-m 



The response of an extended system to an electric 
field is an intrinsic property of central importance to 
the understanding and characterization of bulk materi- 
als. First, it depends very sensitively on an accurate 
description of the correlations between interacting elec- 
trons. Second, it allows one to predict observable quan- 
tities such as the infrared or (non-resonant) Raman and 
hyper-Raman spectra, establishing a direct link between 
macroscopic experimental observations and microscopic 
calculations. Density-functional approaches to calculate 
dielectric susceptibilities in periodic systems were intro- 
duced more than two decades ago However, 
the resulting agreement with experiments is less strik- 
ing than for the case of structural or vibrational prop- 
erties, even though the response to a static electric field 
is strictly a ground-state property 0. These discrep- 
ancies have been related to the intrinsic dependence of 
the exchange and correlation functional on the polariza- 
tion 0; non-local exchange-correlation functionals such 
as the weighted-density approximation (WDA) ^ also 
show some improvements upon the local-density (LDA) 
or generalized-gradient approximations (GGAs) 0. 
While in the case of crystalline semiconductors dielectric 
susceptibilities are usually overestimated by 10 — 30% in 
LDA or GGAs, linear dielectric susceptibilities of con- 
jugated polymers are usually overestimated by a factor 
> 2 — 3, and non- linear susceptibilities by orders of mag- 
nitude (sl llflj . Linear chains of hydrogen dimers are 
another clear example of LDA or GGAs failures, and have 
thus become a stringent test case to assess new develop- 
ments in density-functional theory 

O Hi El 121 Cor- 
related quantum-chemistry methods are able to predict 
accurate linear and non-linear polarizabilities for poly- 
meric chains but their less favourable scaling and 
the need, in the absence of periodic boundary-conditions 
(PBGs), to extrapolate the results to the infinite system 
limit their range of applicability. 

Highly-accurate ground-state electronic-structure cal- 



culations can be performed using continuous quantum 
Monte Garlo (QMC) methods Il5l . which benefit from 
scaling costs as low as TV — [ig , where N is the num- 
ber of electrons. Although diffusion QMC has been suc- 
cessfully applied to address the polarizability of small 
molecules JJJ, the extension to bulk materials is hin- 
dered by the difficulty in treating the response to a lin- 
ear electric field, since the latter is non-periodic and in- 
compatible with PBGs. Only the response of the elec- 
tron gas to a finite wavelength perturbation has been 
investigated In density- functional calculations, the 
response to an electric field has been calculated via its 
long-wavelength limit using supercells of increasing size 
P,|aE^, or, more elegantly, using linear-response the- 
ory [J, 20]. A general approach able to deal with finite 
fields, and thus both linear and non-linear effects, has 
also been recently introduced 0, H^l ■ This method is 
based on the fact that the macroscopic polarization can 
be properly defined even in PBGs |23ll24| : then, instead 
of applying a field in the form a linear potential, a Legen- 
dre transform allows one to introduce an electric enthalpy 
functional whose minimization leads to the correct wave- 
functions. 

In this work we apply these ideas to the case of corre- 
lated electrons, using the many-body formulation of the 
Berry-phase polarization ^25. 26] to deal with the general- 
ization of the electric-enthalpy functional |2l|, |23, 123, |2^ 
to the interacting case jl^. We show how such an ap- 
proach can be applied to the case of diffusion quantum 
Monte Garlo (DMG) calculations, where a local Hermi- 
tian operator is needed to evolve a population of walkers 
representing the ground state. 

Let's consider a system of N interacting electrons in a 
periodic cell of size L (to simplify the notation we de- 
scribe the one-dimensional case, but the extension to 
higher dimensions is straightforward). The normalized 
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many-body wavefunction ^' obeys PBCs 

^'(Xi,..,^^ +L,..,Xn) = 'if{xi,..,Xi,..,XN). (1) 

The polarization of the system can be obtained from the 
single-point Berry phase |25j 
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(2) 
(3) 



where is the size of the cell .3Q] , e is the electron charge, 
G = 2n/L, and X the TV-body operator 



X = Xi + X2 



- XN- 



(4) 



The definition of Eq. |5J) coincides, in the thermody- 
namic limit, with the exact many-body observable j25l |. 
but it also remains well-defined for any finite L. With 
the many-body polarization well defined, the ground- 
state wavefunction of a periodic, extended system in 
the presence of an electric field £ can be obtained 
from the minimum of the generalized electric enthalpy 
I2in2^,27. 28, 29] 



(5) 



where E'^ [5*] is the energy functional for the unperturbed 
Hamiltonian H'^. 

The functional in Eq. |SJ| could be directly minimized 
using a variational Monte Carlo approach (sJI, but it 
can't be applied directly to the more accurate case of 
DMC calculations, since these can only deal with Hamil- 
tonians that are both local and Hermitian. A local Her- 
mitian operator can nevertheless be derived from the 
minimum condition for the electric enthalpy 
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(6) 



with A an appropriate Lagrange multiplier [33 • The ^' 
that minimizes Eq. (0) is also the ground-state for the 
many-body self-consistent Hamiltonian 



H{z) = H° + £^lm— 

2lT z 



iGX 



(7) 



where z depends on 'f through Eq. (|3Jl- For any given L, 
there exists a well-defined interval for £ centered around 
for which the electric enthalpy can be minimized without 
encountering runaway solutions 0, . In addition, all 
physical observables of interest here (e.g. the linear and 
nonlinear dielectric susceptibilities) are derivatives of the 
electric enthalpy with respect to the field and thus remain 
well-defined for every L. 

Due to the self-consistent nature of the operator H{z) 
defined in Eq. the ground state in the presence of an 
electric field must be found through an iterative proce- 
dure. We start from a first value zi for z, e.g. as found in 



DFT-GGA 102.0 123.4 133.5 
DMC 52.2 ± 1.3 55.4 ± 1.2 53.4 ± 1.1 

TABLE I: Linear polarizability per H2 unit for a periodic 
linear chain of H2 dimers, calculated in PBCs for supercells 
containing 10, 16, and 22 H2 units. 



the single-particle calculations or in the many-body trial 
wave function $t; the local Hamiltonian H{zi) is then 
constructed. DMC evolution using this operator leads to 
a new expectation value for z, called Z2, which in turn 
determines a second Hamiltonian H{z2). In the absence 
of stochastic noise, this process could be iterated to con- 
vergence zi ^ Z2 — > 23 ... Zn, to find the fixed 
point of the complex-plane map 



f{Zi) = Zi+l 



(8) 



Since the Monte Carlo procedure introduces a statisti- 
cal error in every estimate of Zi, the map / becomes a 
stochastic function in the complex plane. Now, the lin- 
ear approximation /(z) = az + b is made close to the 
fixed point of f{z) (linearity will be shown numerically 
later). Then, the average over a sequence of {zi\ pro- 
vides the estimate for the fixed point b/{l — a) of /(z), 
from which the polarization is then obtained via Eq. ^ . 
This descends straightforwardly from the consideration 
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is the fixed point of the map 

/(z), since ^ E^ ^» = ^E./(^»-i) = ^ E» /(^«) = 
Ei ("^-^i + ^) (the second equality is valid in the limit 
N 00). Finally, supposing that the error in the com- 
plex plane is isotropic with a quadratic spread a^, the 
{zi} will be distributed around the fixed point with a 
quadratic spread cr^/(l — |ap). 

Most implementations of DMC use importance sam- 
pling: during diffusion in imaginary time, the walkers 
at a time t are distributed as <I't(R)^'(R.), where 5* is 
the correct ground-state wavefunction and <I>t is the trial 
wavefunction (<I>t is chosen here to be the product of a 
Jastrow factor and a Slater determinant). Since the op- 
erator e^'^'^ does not commute with the Hamiltonian, we 
use a forward- walking strategy [s^ [s^ to sample the ex- 
pectation value z. In a branching DMC approach the 
estimate of 4'(Rj.t)/$T(R-j,t) for any walker j at a time 
t and position Rj t is given by the number of its descen- 
dants after a projection time At. The expectation value 
z of e**^^ then becomes 
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(9) 



where is the number of walkers at timestep t, and 
-^j,T-At corresponds to the ascendant configuration of 
the walker j at an interval At back in imaginary time. 

We implemented this approach in the QMC Casino 
package [33; the single-particle orbitals for the Slater 
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Method 


a (a.u.) 


Scaling cost 


PBE-GGA 


144.5 




CCD 


47.6 




CCSD 


48.0 


N« 


CCSD(T) 


50.6 




MP2 


58.0 




MPS 


54.3 


N« 


MP4 


53.6 




DMC 


53.4± 1.1 





TABLE II: Linear polarizability per H2 unit for a periodic 
linear chain of H2 dimers calculated using several methods, 
together with their associated scaling costs. The quantum 
chemistry results are from Ref. [i^ . 



determinant of the trial wavefunction $t were obtained 
from Quantum ESPRESSO jsi], using consistently the 
same finite electric field £ of the QMC calculation |39j . 

To validate our approach, we first calculate the polar- 
izability of an isolated hydrogen atom, using a large su- 
percell and periodic-boundary conditions. For this case, 
a homogeneous finite field can be explicitly applied us- 
ing a saw-tooth potential (i.e. a potential that is both 
periodic and piecewise linear), and the results compared 
with those obtained with our method. We use the same 
GGA trial wavefunctions for both tests, as calculated un- 
der an applied electric field of 0.005 a.u. using a local H 
pseudopotential, 50 Ry cutoff in the plane- wave basis set, 
and a cell of 20 x 20 x 200 a. u. (the longest direction 
coincides with the applied electric field, and a large su- 
percell has been used to eliminate finite-size errors |4ojV 
In all DMC calculations we choose a timestep of 0.02 
a.u., assuring an acceptance ratio greater than 99.7%. 
Then, the projection time At for the forward-walking is 
determined. We show in Fig.ma) the dependence of the 
expectation value for the polarization of the atom in the 
absence of an applied field, as a function of the projec- 
tion time At in Eq. lO, showing convergence to a.u. 
for At > 600 timesteps. With these parameters, we start 
our iterative sequence of calculations. The first value for 
z is taken from the GGA results, and then after each 
DMC run the local operator H{z) is updated using the 
current expectation value for z. We report in Fig. ^b) 
the atomic polarizability calculated from each one of the 
Zi along the sequence of DMC runs. The average of the 
Zi over the last 13 runs gives us a fixed point estimate 
(0.9995072 ±2.* IQ-'^ + 7.05 * lO^S ± 4. * IQ-^i) for 
/(z), and a corresponding polarizability for the H atom 
of 4.49 ±0.03 a.u., in excellent agreement with the result 
obtained for the saw-tooth potential (4.52 ± 0.05 a.u.) 
and the exact value of 4.5 a.u. 0- 

Having validated our scheme in an isolated case, 
we proceed to apply it to a genuine periodic system, 
namely a linear chain of hydrogen dimers, where stan- 
dard density-functional theory performs very poorly [lC| | . 
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FIG. 1: (a) Polarization of an isolated H atom in the absence 
of an electric field as a function of the forward walking param- 
eter At, using a trial wavefunction expressed by a PBE-GGA 
orbital calculated in the presence of an electric field of 0.005 
a.u. . The dashed line is a guide to the eye. (b) Polarizability 
of an isolated H atom, for each of the self-consistent DMC 
runs in the iterative sequence (see text). The dashed line is a 
guide to the eye. 
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FIG. 2: (a) Linear polarizability per H2 unit for a periodic 
linear chain of 22 H2 dimers in PBCs, for each of the self- 
consistent DMC runs in the iterative sequence (see text). The 
dashed line is a guide to the eye. (b) Last ten Zi values (thin), 
and estimate of the fixed point of f{z) (bold). 



We adopt one of the standard geometries reported in 
the literature 0|, where the hydrogen dimers have 
a H-H bond length of 2.0 a.u. and are separated by 2.5 
a.u. . For this case, reliable reference results for the 
polarizability per dimer unit a are available from cor- 
related quantum-chemistry approaches (a = 53.6 a.u. 
in MP4 and a = 50.6 a.u. in CCSD(T)) [i^; these 
numbers should be contrasted with the PBE-GGA result 
a = 144.5 a.u. ji^l. To calculate the DMC suscepti- 
bility, we use a trial wavefunction composed of a Slater 
determinant of PBE-GGA single-particle orbitals (con- 
sistently calculated with an electric field of 0.001 a.u.) 
times a Jastrow term containing the electron-electron and 
electron- nucleus terms 01 • We consider three supercells 
containing 10, 16 and 22 H2 units, and the same elec- 
tric field of 0.001 a.u. . The polarizabilities per unit 
dimer are found by finite differences (note that the po- 
larization at zero field is zero by symmetry). In Fig. 
I3a), we show the polarizability calculated at each DMC 
run in the 22-unit supercell. After the second itera- 
tion, the instantaneous polarizability oscillates roughly 
around the final value; again, we stress that it is the av- 
erage of {zi} that delivers the estimate of the fixed point 
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and thus the final polarization. In Fig. Efb), we show 
the distribution of the last ten Zi, and the fixed point 
(0.8029±4.*10-'*+6.0*10~2i±i.*10-3i), obtained from 
their average. In Table U we illustrate the convergence 
of the calculated polarizability with respect to the super- 
cell size ^3 : The faster convergence with respect to the 
PBE-GGA case originates in the the stronger localization 
of 'J in DMC. Indeed, the localization spread (defined 
as -(i2/iV47r2)ln|zp ji^) decreases from 4.32 a.u. (i.e. 
Bohr^) for PBE-GGA to 2.44 ±0.01 a.u. in DMC for the 
22 units supercell. Finally, we compare in Table ^1 our 
final result a — 53.4 ± 1.1 a.u. with those obtained using 
other approaches: we find excellent agreement with the 
most accurate results available using MP3 and MP4. We 
note in passing that the correlated quantum chemistry 
polarizabilities have been extrapolated from calculations 
performed for finite chains |43, with the extrapolative 
correction amounting to 20% of the final results. 

In conclusion, we introduced and validated a novel ap- 
proach to study periodic systems in the presence of fi- 
nite electric-fields with highly-accurate, and in principle 
linear-scaling, diffusion QMC calculations. In particular, 
we presented the first QMC evaluation of the dielectric 
susceptibility of an extended system. For the paradig- 
matic test case of hydrogen chains, we obtained excellent 
agreement with the best quantum-chemistry results. Our 
approach opens the possibility to evaluate linear and non- 
linear dielectric properties, effective charges, and Raman 
and infrared coupling tensors using QMC calculations. 
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